home *** CD-ROM | disk | FTP | other *** search
/ MacHack 1997 / MacHack 1997.toast / Hacks / Hacks ’95 / Venus / project_3D.cc < prev    next >
Text File  |  1995-06-23  |  18KB  |  539 lines
  1. // This may look like C code, but it is really -*- C++ -*-
  2. /*
  3.  ************************************************************************
  4.  *
  5.  *                           Rendering of a 3D function
  6.  *
  7.  * Suppose we have a function z=f(x,y) specified as an image, where
  8.  * pixel intensity z at point (x,y) is the value of the function f(x,y).
  9.  * We assume the map (function) is periodic, i.e., f(x+n*Dx,y+m*Dy)=f(x,y)
  10.  * for each integer n and m; Dx and Dy are dimensions of the map. We assume
  11.  * Dx = Dy = D. In the following, we assume that x, y are always within
  12.  * [0,D-1].
  13.  *
  14.  * The coordinate system (x,y,z) is chosen so that x increases along
  15.  * the horizontal lines, z is the elevation, and y increases "in depth".
  16.  *
  17.  * We would use a perspective projection with the center (x0, y0, z0)
  18.  * and the projection plane y=y0+eye_focus. So, the viewer looks from the
  19.  * point (x0,y0,z0) straight "in depth".
  20.  * When performing projection, we project the points of the map with y within 
  21.  *         [y0+eye_focus,y0+sight_depth]
  22.  * In doing the projections, we scan the 3D map along the line with y=const
  23.  * from y=y0+sight_depth to y=y0+eye_focus, that is, from farther to closer.
  24.  *
  25.  * Thus, we need to project a volume {x in [0,D), y in [y0+eye_focus,y0+sight_depth],
  26.  * z in [0,Z)} into the view plane {u in (0,Du), v in (0,Dv)}
  27.  * View plane is associated with the eye focal plane (and fixed at the observer,
  28.  * that is, us). We use a perspective transformation: draw a line from a point
  29.  * (x1,y1,z1) of the volume to a view point (x0,y0,z0), and see where the line
  30.  * intersects the view plane y=y0+eye_focus.
  31.  * The equation for the line connecting (x1,y1,z1) with (x0,y0,z0) is
  32.  *    (x-x0)/(x1-x0) = (y-y0)/(y1-y0) = (z-z0)/(z1-z0)
  33.  * The line intersects the plane y=y0+eye_focus at a point
  34.  *    x = x0 + (x1-x0)*factor
  35.  *  z = z0 + (z1-z0)*factor
  36.  *  with factor = eye_focus/(y1-y0)
  37.  * To finish the translation into the view port coordinates, we need to know
  38.  * the projection of our eye: if we look straight into the horizon, that is,
  39.  * y1=infinity, we see a point (x0,y0) at the view plane. We want this point
  40.  * to be suitably located in the view plane, that is, having the view plane
  41.  * coordinates (u=Du/2, v=horizon). It gives us necessary transformation formulas
  42.  * from a point (x,y,z) into the point (u,v) of the view plane
  43.  *    u = Du/2 + (x-x0)*factor
  44.  *  v = horizon + (z-z0)*factor
  45.  *    where
  46.  *        factor = eye_focus/(y-y0)
  47.  * Note that y >= y0+eye_focus always, therefore y>y0. That is, we can't see
  48.  * (at least distinctly) objects very close to our eyes, let alone objects
  49.  * behind our eyes. So, for all y we consider, 0<factor<=1
  50.  *
  51.  * If we scan the map along the line y=const, then 'factor' remains constant
  52.  * during the scan. Of course, we can increment x by one and then find
  53.  * the corresponding point (u,v) on the view plane and draw it.
  54.  * However, since the view plane is what we are going to display, it's
  55.  * better to increment u by one and then work it out back to x, find
  56.  * z value at (x,y) and compute v.
  57.  *
  58.  * So, the algorithm to process the scanline y=const is as follows:
  59.  * 1. fix y=const within [y0+eye_focus,y0+sight_depth] and compute
  60.  *    factor=eye_focus/(y-y0)
  61.  * 2. compute xref = -D/2/factor + x0
  62.  * 3. for u=0 to D-1
  63.  * 4.     compute x=(round)xref
  64.  * 5.     find z=f(x,y) from the map
  65.  * 6.      compute v = factor*(z-z0) + horizon
  66.  * 7.     xref += 1/factor
  67.  *     endfor
  68.  *
  69.  * Note, that the 'factor' should be represented as a fraction, since 0<factor<=1
  70.  * Division and multiplications by the factor are fixed-point arithmetic
  71.  * operations; xref is also a fixed-point number, the others are integers.
  72.  * We also assume that Du (width of the view plane) is equal to D (width of
  73.  * the map), though it isn't so critical.
  74.  *
  75.  * Since we scan from larger y's to smaller y's (that is, from farther ahead
  76.  * to near, that is, towards the observer), then if a (u,v) point generated by
  77.  * the algorithm should obscure (u,v) point generated at an earlier pass
  78.  * (with larger y), it would just overwrite the old point. However, if we just
  79.  * draw (u,v) dots we compute we get a dotty picture. Note, that actually we see
  80.  * a segment of line (x,y,z)-(x1,y-dy,z1). If we just connect two points generated
  81.  * in two consecutive scans, (u,v1) and (u,vold), it wouldn't be always
  82.  * right, because the original line could've been obscured. Suppose that
  83.  * z0 (elevation of the observation point) is high enough (comparable with
  84.  * the hight of the tallest peaks). Then if v(u,y+dy) > v(u,y) for some
  85.  * u and dy>0, then we have a descending slope, and the line is visible.
  86.  * Otherwise the line is on ascending slope, and it couldn't be visible
  87.  * because it will be obscured by the descending slope (which should be
  88.  * closer to us). Note the line from (u,c1) to (u,vold) is a vertical one,
  89.  * and therefore, is very easy to draw (and clip, if necessary).
  90.  *
  91.  * When we create (u,v) map, we assign to the point (u,v) the intensity
  92.  * (color) of the f(x,y) point of the original map. We can use Gouraud
  93.  * interpolation in drawing the line.
  94.  * Generalization: draw boxes (quadraluzation), than we can increase the
  95.  * scanning step in u (or in y) in the algorithm above.
  96.  *
  97.  * For this particular program, we also draw a part of the map as it is
  98.  * (that is, a 2D picture) above the horizon: just to make clouds.
  99.  *
  100.  * Inspiration:
  101.  * Tim Clarke, tjc1005@hermes.cam.ac.uk
  102.  * Note, there are some typos in that post. And my notation is completely different
  103.  * from his.
  104.  *
  105.  ************************************************************************
  106.  */
  107.  
  108. #include "std.h"
  109. #include "image.h"
  110. #include "window.h"
  111.  
  112. #define Debug 0
  113. const int eye_focus = 70;
  114. const int sight_depth = 170;
  115.  
  116. /*
  117.  *----------------------------------------------------------------------
  118.  *                An instance of a generic screen window to display
  119.  *                                 a "3D" image
  120.  */
  121.  
  122. class ProjectionWindow : public OffScreenWindow
  123. {
  124.     const IMAGE& map;                    // that specifies z=f(x,y)
  125.  
  126.     int x0, y0, z0;                    // Point of view
  127.     int hor;                            // Elevation of the horizon on the view plane
  128.     Boolean do_animation;
  129.     
  130.     void project(void);                    // Draw a projection in an offscreen
  131.                                         // world
  132.                                         // Draw clouds above the horizon
  133.     void draw_clouds(const int horizon, const int starting_y);
  134.     
  135.                                 // A class to handle one scanline of the view plane
  136.     struct OneViewScanline
  137.     {
  138.         const int width;
  139.         short * vs;
  140.         char  * color;
  141.         
  142.         OneViewScanline(const int _width);
  143.         ~OneViewScanline(void);    
  144.                                 // Put all zeros (for the final scanline)
  145.         void clear(void);
  146.         
  147.     };
  148.     
  149.     OneViewScanline scanline1, scanline2;
  150.     void draw_scanline(OneViewScanline& scanline, const int y);
  151.                                 // Connect two scanlines, sl0 to sl1, drawing
  152.                                 // vertical lines between corresponding points
  153.                                 // (if the lines are visible)
  154.     void connect_scanlines(OneViewScanline& sl0, OneViewScanline& sl1);
  155.  
  156.                                         // redefining some event handlers
  157.     virtual Boolean handle_key_down(const EventRecord& the_event);
  158.     virtual Boolean handle_null_event(const long event_time);
  159.     
  160. public:
  161.     ProjectionWindow(const IMAGE& image, const char * title);
  162.     ~ProjectionWindow(void)  {}
  163. };
  164.  
  165.  
  166.                                 // Creating a window that would have our picture
  167.                                 // displayed.
  168. ProjectionWindow::ProjectionWindow(const IMAGE& image, const char * title)
  169.     : OffScreenWindow(ScreenRect(rowcol(256,image.q_ncols())),title,256),
  170.       map(image),
  171.       scanline1(image.q_ncols()),
  172.       scanline2(image.q_ncols()),
  173.       do_animation(FALSE)
  174. {
  175.     x0 = map.q_ncols()/2;            // Pick up the initial observation point
  176.     y0=0;                            // somewhere in the middle
  177.     z0 = 200;
  178.     hor = 150;
  179.     
  180.     project();
  181.     do_animation = TRUE;
  182. }
  183.  
  184.  
  185.                                     // Handles key_down & auto_key events. Return FALSE
  186.                                     // if the window is to be closed down
  187.                                     // The key handled move the observer and/or
  188.                                     // change his direction of view (horizon)
  189.                                     // To the observer (that is, us) it looks like
  190.                                     // the entire scene moves.
  191. Boolean ProjectionWindow::handle_key_down(const EventRecord& the_event)
  192. {
  193.     //message("char pressed %ld",the_event.message & charCodeMask);
  194.     switch(the_event.message & charCodeMask)
  195.     {
  196.         case 11:                // PgUp
  197.              if( the_event.modifiers & optionKey )
  198.                if( hor < 250 )
  199.                  hor++;
  200.                else if( z0 < 250 )
  201.                  z0++;
  202.              break;
  203.              
  204.         case 12:                // PgDn
  205.              if( the_event.modifiers & optionKey )
  206.                if( hor > 5 )
  207.                  hor--;
  208.                else if( z0 > 0 )
  209.                  z0--;
  210.              break;
  211.              
  212.         
  213.         case 30:                // Arrow Up
  214.              if( (y0 +=2) >= map.q_nrows() )
  215.                y0 = 0;
  216.              break;
  217.              
  218.         case 31:                // Arrow Down
  219.              if( (y0 -=2) < 0 )
  220.                y0 = map.q_nrows()-1;
  221.              break;
  222.  
  223.         case 28:                // Arrow Left
  224.              if( (x0 -=4) < 0 )
  225.                x0 = map.q_ncols()-1;
  226.              break;
  227.  
  228.         case 29:                // Arrow Righ
  229.              if( (x0 +=4) >= map.q_ncols() )
  230.                x0 = 0;
  231.              break;
  232.              
  233.         default:
  234.              return FALSE;                // Other keys kill the application
  235.              
  236.     }
  237.     project();
  238.     refresh();
  239.     return TRUE;
  240. }
  241.  
  242.  
  243.                                     // Handles null events, when nothing happens for some
  244.                                     // time. Return FALSE when it's time to die
  245.                                     // This function changes the view point (that
  246.                                     // is, moves the scene) with the time
  247.                                     // (or with the mouse moves), and makes animation
  248.                                     // When the mouse button is pressed, the
  249.                                     // scene moves with the mouse; otherwise, it
  250.                                     // moves by itself.
  251. Boolean ProjectionWindow::handle_null_event(const long event_time)
  252. {
  253.   static long int prev_tick_count = 0;
  254.   static Point prev_mouse_loc;
  255.   
  256.   if( prev_tick_count == 0 )
  257.   {
  258.     prev_tick_count = TickCount();
  259.     GetMouse(&prev_mouse_loc);
  260.     return TRUE;
  261.   }
  262.   
  263.   if( !do_animation )
  264.     return TRUE;
  265.  
  266.   if( StillDown() )                    // If mouse button is kept pressed, the user
  267.   {                                    // wants to lead the way. Move the picture
  268.       Point new_point;                // as he moves the mouse
  269.       GetMouse(&new_point);
  270.       x0 += (new_point.h - prev_mouse_loc.h)*2;
  271.       y0 -= (new_point.v - prev_mouse_loc.v)*2;
  272.       prev_mouse_loc = new_point;
  273.   }
  274.   else
  275.   {
  276. //    prev_tick_count = TickCount();
  277.  
  278.     x0 += 4;                                // Move the scene - do animation
  279.     y0 += 4;
  280.   }
  281.   
  282.   if( x0 >= map.q_ncols())                    // Do some clipping
  283.     x0 = 0;
  284.   else if( x0 < 0 )
  285.     x0 = map.q_ncols()-1;
  286.   if( y0 >= map.q_nrows())
  287.     y0 = 0;
  288.   else if( y0 < 0 )
  289.     y0 = map.q_nrows()-1;
  290.  
  291.   project();
  292.   refresh();
  293.     
  294.   return TRUE;                            // Keep going
  295. }
  296.  
  297.  
  298.  
  299.                                 // Allocate data for one scanline of the view window
  300. ProjectionWindow::OneViewScanline::OneViewScanline(const int _width)
  301.     : width(_width)
  302. {
  303.                                     // Allocate memory in one chunk for both arrays
  304.     vs = (short *)malloc(width*(sizeof(vs[0])+sizeof(color[0])));
  305.     assert( vs != 0 );
  306.     color = (char *)vs + width*sizeof(vs[0]);
  307. }
  308.  
  309.                                 // Dispose of the arrays
  310. ProjectionWindow::OneViewScanline::~OneViewScanline(void)
  311. {
  312.     assert( vs != 0 );
  313.     delete vs;                        // color would be disposed of, too
  314. }
  315.  
  316.                                 // Put all zeros (for the final scanline)
  317. void ProjectionWindow::OneViewScanline::clear(void)
  318. {
  319.     memset((void *)vs,0,width*sizeof(vs[0]));
  320.     memset((void *)color,0,width*sizeof(color[0]));
  321. }
  322.  
  323.  
  324.                                 // Draw a scanline, i.e. line of const y for
  325.                                 // all u's within 0..width-1
  326.                                 // see formulas above
  327.                                 // For factor, invfactor, xref we use a fixed-
  328.                                 // point arithmetic with 8-bit implied fraction
  329. void ProjectionWindow::draw_scanline(OneViewScanline& scanline, const int y)
  330. {
  331.     int factor = (eye_focus<<8)/(y-y0);
  332.     int invfactor = ((long)(y-y0)<<8)/eye_focus;
  333.                                     // We know that xref = -D/2/factor + x0
  334.                                     // But we want to keep xref positive, moreover,
  335.                                     // within [0,D-1]. Since the map is periodic
  336.                                     // with period D, then we can write
  337.                                     // xref = floor(1/2/factor)*D - (D/2)*(1/factor) + x0
  338.                                     // or,
  339.                                     // xref = D( floor(1/2/factor) - (1/2/factor) ) + x0
  340.                                     //      = -D*frac(1/2/factor) + x0
  341.                                     // In fixed point 8-bit fraction arithmetic, it's elementary
  342.                                     // to separate integer and fractional part of
  343.                                     // the invfactor/2
  344.                                     // It's easy to see that xref we got is within the desired
  345.                                     // interval, give or take one D.
  346.     int xref = (x0<<8) - scanline.width * ( (unsigned char)(invfactor >> 1) );
  347.                                     // Make sure xref is within [0,width) in
  348.                                     // the fixed-point 8-bit fraction arithmetic
  349.     const int width_fp = scanline.width<<8;
  350.     if( xref >= width_fp )
  351.         xref -= width_fp;
  352.     if( xref < 0 )
  353.         xref += width_fp;
  354.      
  355.     
  356.                                         // Clip y to the map, if map.q_nrows()
  357.                                         // is an exact power of two, clipping
  358.                                         // is just &
  359.     int y_map = y & (map.q_nrows()-1);
  360. #if 0
  361.     int y_map = y;
  362.     while( y_map >= map.q_nrows() )
  363.         y_map -= map.q_nrows();
  364. #endif
  365.  
  366. #if Debug
  367.     message("draw scanline y=%d, factor %d, xref %d",y,factor,xref);
  368. #endif
  369.     register short * vp;
  370.     register char * cp;
  371.     for(vp=scanline.vs, cp = scanline.color; vp<&scanline.vs[scanline.width]; )
  372.     {
  373.         int x = xref>>8;
  374.         int z = map(y_map,x);
  375.         *vp++ = ((factor*(z-z0))>>8) + hor;
  376.         *cp++ = z;
  377. #if Debug
  378.         if( vp - scanline.vs < 5 )
  379.           message("map(%d,%d) is %d, projected to v=%d",y_map,x,z,*(vp-1));
  380. #endif
  381.         if( (xref += invfactor) >= width_fp )
  382.           xref -= width_fp;
  383.     }
  384. }
  385.  
  386.  
  387.                                     // Draw a projection in an offscreen
  388.                                     // graf world
  389. void ProjectionWindow::project(void)
  390. {
  391. //    SetOffscreenWorld offscreen_world(*this);
  392.  
  393. //    ScreenRect sub_horizon(q_bounds());
  394. //    sub_horizon.top +=40; //= sub_horizon.bottom-z0;
  395. //    sub_horizon.print("sub horizon");
  396. //    EraseRect(sub_horizon);
  397.  
  398.     draw_clouds(hor,eye_focus);
  399.  
  400.     OneViewScanline * curr_scanline = &scanline1,
  401.                     * prev_scanline = &scanline2;
  402.                                                 // note curr_scanline ^ pointer_diff
  403.                                                 // gives prev_scanline, and
  404.                                                 // vice versa
  405.     int pointer_diff = (long int)curr_scanline ^ (long int)prev_scanline;
  406.     
  407.     draw_scanline(*prev_scanline,y0+sight_depth);
  408.  
  409.     register int y;
  410.     for(y=y0+sight_depth-1; y >= y0 + eye_focus; y--)
  411.     {
  412. #if Debug
  413.         if((y0+sight_depth-1)-y > 3 )
  414.           exit(0);
  415. #endif
  416.         draw_scanline(*curr_scanline,y);
  417.         
  418.         connect_scanlines(*prev_scanline,*curr_scanline);
  419.         
  420.                                             // exchange pointers
  421. //        curr_scanline = (OneViewScanline*)((long int)curr_scanline ^ pointer_diff);
  422. //        prev_scanline = (OneViewScanline*)((long int)prev_scanline ^ pointer_diff);
  423.         (long int&)prev_scanline ^= pointer_diff;
  424.         (long int&)curr_scanline ^= pointer_diff;
  425.     }
  426. }
  427.  
  428.                                 // Connect two scanlines, sl0 to sl1, drawing
  429.                                 // vertical lines between the corresponding points
  430.                                 // (if the lines are visible)
  431. void ProjectionWindow::connect_scanlines(
  432.     OneViewScanline& sl0, OneViewScanline& sl1)
  433. {
  434.     PixMapHandle pixmap = get_pixmap();
  435.  
  436.     assert( LockPixels(pixmap) );
  437.     
  438.     char * pixp = GetPixBaseAddr(pixmap);
  439.     register int u;
  440.     for(u=0; u<width(); u++)                    // Draw a vertical line (u,v)-(u,vold)
  441.     {                                            // (providing it's visible)
  442.       int vold = sl0.vs[u];
  443.       int v    = sl1.vs[u];
  444.       if( v > vold )
  445.         continue;                    // The line is on the ascending slope and *will* be obscured later
  446.                                     // Keep in mind that v <= vold now
  447.       if( v >= height() || vold <= 0 )
  448.         continue;                    // means both v and vold are out-of-picture
  449.  
  450.       if( v < 0 )
  451.         v = 0;
  452.                                     // Now, 0<=v <= vold
  453.       int z0    = (unsigned char)sl0.color[u]<<4;    // For the purpose of color interpolation
  454.       int z     = (unsigned char)sl1.color[u]<<4;    // we use fixed-point arithmetic with 4 bit
  455.       int stretch = v == vold ? 0 : (z0-z)/(vold-v);// implicit fraction
  456. #if Debug
  457.     if(u<4)
  458.       message("u=%ld drawing from v=%ld to vold=%ld at zs %lx-%lx",u,v,vold,z>>4,z0>>4);
  459. #endif
  460.                                       // draw a line from v to vold
  461.                                       // Note that v on mac goes from top to bottom,
  462.                                       // that is, we need a flip
  463.       char * pixp_beg = pixp + u + (height()-1-v) * bytes_per_row();
  464. //      *pixp_beg = 126;
  465.       for(; v <= (vold > height()-1 ? height()-1 : vold); v++,pixp_beg -= bytes_per_row())
  466.           *pixp_beg = z>>4,
  467.           z += stretch;
  468.     }
  469.  
  470.     UnlockPixels(pixmap);
  471. }
  472.  
  473.                                 // Draw clouds above the horizon line
  474.                                 // using the same map
  475.                                 // Simply speaking, just copy a rectangular
  476.                                 // segment of the map with y from 0 to
  477.                                 // horizon into the offscreen world
  478.                                 // Below the horizon, fill everything (I mean,
  479.                                 // the offscreen world) with the
  480.                                 // background "pixel"
  481. void ProjectionWindow::draw_clouds(const int horizon, const int starting_y)
  482. {
  483.     const unsigned char background_pixel = 124;
  484.     
  485.     assert( horizon >= 0 && horizon < height() );
  486.     
  487.     PixMapHandle pixmap = get_pixmap();
  488.     
  489.     assert( LockPixels(pixmap) );
  490.     
  491.     char * pixp = GetPixBaseAddr(pixmap);
  492.     char * pixp_stop = pixp + (height()-horizon)*bytes_per_row();
  493.     register int i;                                // "Draw" the image
  494.     for(i=starting_y; pixp<pixp_stop; i++, pixp += bytes_per_row()-width())
  495.     {                                            // Note, bytes_per_row is generally GREATER
  496.       register int j;                            // than width, and is always a multiple of 4
  497.       if( i >= map.q_nrows() )
  498.         i = 0;                                    // make wrap-around
  499.       for(j=0; j<width(); j++)
  500.            *pixp++ = map(i,j);
  501.     }
  502.  
  503.     pixp_stop += horizon*bytes_per_row();        // Continue till the end of pixmap
  504.     for(; pixp<pixp_stop; pixp += bytes_per_row()-width())
  505.     {                                            // Note, bytes_per_row is generally GREATER
  506.       register int j;                            // than width, and is always a multiple of 4
  507.       for(j=0; j<width(); j++)
  508.            *pixp++ = background_pixel;
  509.     }
  510.  
  511.     UnlockPixels(pixmap);
  512. }
  513.  
  514.  
  515. /*
  516.  *----------------------------------------------------------------------
  517.  *                        Routing display function
  518.  */
  519.  
  520. void project_3D(const IMAGE& image, const char * title)
  521. {
  522.   if( image.q_nrows() & (image.q_nrows()-1) )
  523.     _error("Number of rows, %d isn't an exact power of two",image.q_nrows());
  524.   if( image.q_ncols() & (image.q_ncols()-1) )
  525.     _error("Number of columns, %d isn't an exact power of two",image.q_ncols());
  526. #if 0
  527.     IMAGE image1(image);
  528.     register int i,j;
  529.     int ncols = image.q_ncols();
  530.     for(j=-20; j<20; j++)
  531.       for(i=-5; i<5; i++)
  532.         image1(eye_focus+i,j+ncols/2) = 128 - 5*abs(i+j),
  533.         image1(eye_focus+40+i,j+ncols/2) = 180 - 7*abs(i+j);
  534.     ProjectionWindow map_projection(image1,title);
  535. #else
  536.     ProjectionWindow map_projection(image,title);
  537. #endif
  538.     map_projection.handle();
  539. }